home *** CD-ROM | disk | FTP | other *** search
/ Suzy B Software 2 / Suzy B Software CD-ROM 2 (1994).iso / adult_ed / lectures / lectqbas.txt < prev    next >
Text File  |  1995-05-02  |  4KB  |  200 lines
  1. DEFINT I-M
  2. DEFDBL A-H, O-Z
  3. OPTION BASE 1
  4. DIM SHARED ROT(3, 3), XO(3), XN(3), PI2
  5. DECLARE SUB R3 (TH)
  6. DECLARE SUB R2 (TH)
  7. DECLARE SUB R1 (TH)
  8. DECLARE SUB M1 (ROT)
  9. DECLARE SUB M2 (ROT)
  10. DECLARE SUB M3 (ROT)
  11. DECLARE SUB MATMULT (XO, XN, ROT)
  12. DECLARE SUB XHAT (XLON, XLAT, XO)
  13. DECLARE SUB CLR (XO, XN)
  14. DECLARE FUNCTION PLG (X, Y)
  15. PI = ATN(1#) * 4#
  16. DEGRAD = PI / 180!
  17. PI2 = PI / 2#
  18. CLS
  19. '
  20. ' Input angles:
  21. '    PHI = site latitude in DEGREES
  22. '    DEC = declination of star in DEGREES
  23. '    RA = right ascension of star in DEGREES (e.g., 12h = 180 deg)
  24. '    ST = local sidereal time (actual or mean depending on which
  25. '         equinox has been chosen for reference) in DEGREES
  26. '
  27. PHI = 45#: PHIR = PHI * DEGRAD
  28. DEC = 0#: DECR = DEC * DEGRAD
  29. RA = 180#: RAR = RA * DEGRAD
  30. ST = 0#: STR = ST * DEGRAD
  31. '
  32. ' First rotate the RA and DEC (Q system) coordinates into HA and DEC
  33. ' (A system) by the formula X(HA,DEC)=M2*R3(ST)*X(RA,DEC)
  34. '
  35. CALL XHAT(RAR, PI2 - DECR, XO)
  36. CALL R3(STR)
  37. CALL MATMULT(XO, XN, ROT)
  38. CALL CLR(XO, XN)
  39. CALL M2(ROT)
  40. CALL MATMULT(XO, XN, ROT)
  41. CALL CLR(XO, XN)
  42. '
  43. ' Vector XO now contains the components of the star position in the A system
  44. '
  45. ' Now we convert from the A system to the H (horizon) system by the formula
  46. '    X(AZ,ALT)=R3(180)*R2(90-PHI)*X(HA,DEC)
  47. '
  48. ANG = PI2 - PHIR
  49. CALL R2(ANG)
  50. CALL MATMULT(XO, XN, ROT)
  51. CALL CLR(XO, XN)
  52. CALL R3(PI)
  53. CALL MATMULT(XO, XN, ROT)
  54. CALL CLR(XO, XN)
  55. '
  56. ' Vector XO contains the components of the star position in the H system
  57. '
  58. ' Now we need to get the altitude and azimuth angles
  59. '
  60. X = XO(1)
  61. Y = XO(2)
  62. IF (X <> 0# OR Y <> 0#) THEN
  63.   ARG = XO(3) / (SQR(XO(1) * XO(1) + XO(2) * XO(2)))
  64.   ALTR = ATN(ARG)
  65.   ALT = ALTR / DEGRAD
  66.   PRINT "Altitude angle is "; : PRINT USING "###.###"; ALT; : PRINT "degrees."
  67. ELSE
  68.   IF (XO(3) > 0) THEN
  69.     PRINT "Object is at the zenith."
  70.   ELSE
  71.     PRINT "Object is at the nadir."
  72.   END IF
  73. END IF
  74. AZ = PLG(X, Y) / DEGRAD
  75. IF (AZ <= 360#) THEN
  76.   PRINT "Azimuth angle is "; : PRINT USING "###.###"; AZ; : PRINT "degrees."
  77. ELSE
  78.   PRINT "Azimuth is undefined."
  79. END IF
  80. SUB CLR (XO, XN)
  81. '
  82. ' This subroutine replaces vector XO with vector XN
  83. '
  84. FOR I = 1 TO 3
  85.     XO(I) = XN(I)
  86.     TEST = ABS(XO(I))
  87.     IF (TEST < .000000000000001#) THEN
  88.        XO(I) = 0#
  89.     END IF
  90. NEXT I
  91. END SUB
  92. SUB M1 (ROT)
  93. ROT(1, 1) = -1#
  94. ROT(1, 2) = 0#
  95. ROT(1, 3) = 0#
  96. ROT(2, 1) = 0#
  97. ROT(2, 2) = 1#
  98. ROT(2, 3) = 0#
  99. ROT(3, 1) = 0#
  100. ROT(3, 2) = 0#
  101. ROT(3, 3) = 1#
  102. END SUB
  103. SUB M2 (ROT)
  104. ROT(1, 1) = 1#
  105. ROT(1, 2) = 0#
  106. ROT(1, 3) = 0#
  107. ROT(2, 1) = 0#
  108. ROT(2, 2) = -1#
  109. ROT(2, 3) = 0#
  110. ROT(3, 1) = 0#
  111. ROT(3, 2) = 0#
  112. ROT(3, 3) = 1#
  113. END SUB
  114. SUB M3 (ROT)
  115. ROT(1, 1) = 1#
  116. ROT(1, 2) = 0#
  117. ROT(1, 3) = 0#
  118. ROT(2, 1) = 0#
  119. ROT(2, 2) = 1#
  120. ROT(2, 3) = 0#
  121. ROT(3, 1) = 0#
  122. ROT(3, 2) = 0#
  123. ROT(3, 3) = -1#
  124. END SUB
  125. SUB MATMULT (XO, XN, ROT)
  126. '
  127. ' This subroutine computes the product of a vector and a rotation matrix
  128. '
  129. XN(1) = ROT(1, 1) * XO(1) + ROT(1, 2) * XO(2) + ROT(1, 3) * XO(3)
  130. XN(2) = ROT(2, 1) * XO(1) + ROT(2, 2) * XO(2) + ROT(2, 3) * XO(3)
  131. XN(3) = ROT(3, 1) * XO(1) + ROT(3, 2) * XO(2) + ROT(3, 3) * XO(3)
  132. END SUB
  133. FUNCTION PLG (X, Y)
  134.   SGNAX = SGN(ABS(X))
  135.   SGNAY = SGN(ABS(Y))
  136.   SGNX = SGN(X)
  137.   SGNY = SGN(Y)
  138.   C = PI2 * (2# - (1# + SGNX) * (SGNY - SGNAY + 1#))
  139.   IF (X <> 0#) THEN
  140.     PLG = SGNAX * ATN(Y / X) + C
  141.   ELSE
  142.     IF (Y <> 0) THEN
  143.        PLG = C
  144.     ELSE
  145. '      PLG not defined for X and Y both zero.
  146.        PLG = 9999#
  147.     END IF
  148.   END IF
  149. END FUNCTION
  150. SUB R1 (TH)
  151. 'DIM ROT(3, 3)
  152. CTH = COS(TH)
  153. ST = SIN(TH)
  154. ROT(1, 1) = 1!
  155. ROT(1, 2) = 0!
  156. ROT(1, 3) = 0#
  157. ROT(2, 1) = 0#
  158. ROT(2, 2) = CTH
  159. ROT(2, 3) = STH
  160. ROT(3, 1) = 0#
  161. ROT(3, 2) = -STH
  162. ROT(3, 3) = CTH
  163. END SUB
  164. SUB R2 (TH)
  165. 'DIM ROT(3, 3)
  166. CTH = COS(TH)
  167. STH = SIN(TH)
  168. ROT(1, 1) = CTH
  169. ROT(1, 2) = 0#
  170. ROT(1, 3) = -STH
  171. ROT(2, 1) = 0#
  172. ROT(2, 2) = 1#
  173. ROT(2, 3) = 0#
  174. ROT(3, 1) = STH
  175. ROT(3, 2) = 0#
  176. ROT(3, 3) = CTH
  177. END SUB
  178. SUB R3 (TH)
  179. CTH = COS(TH)
  180. STH = SIN(TH)
  181. ROT(1, 1) = CTH
  182. ROT(1, 2) = STH
  183. ROT(1, 3) = 0#
  184. ROT(2, 1) = -STH
  185. ROT(2, 2) = CTH
  186. ROT(2, 3) = 0#
  187. ROT(3, 1) = 0#
  188. ROT(3, 2) = 0#
  189. ROT(3, 3) = 1#
  190. END SUB
  191. SUB XHAT (XLON, XLAT, XO)
  192. '
  193. ' This subroutine computes the unit vector given the longitude and
  194. ' co-latitude angles (in radians)
  195. '
  196. XO(1) = SIN(XLAT) * COS(XLON)
  197. XO(2) = SIN(XLAT) * SIN(XLON)
  198. XO(3) = COS(XLAT)
  199. END SUB
  200.